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Abstract 

We formulate and test a hybrid fluid-Monte Carlo scheme for the treatment of elastic col- 
lisions in gases and plasmas. While our primary focus and demonstrations of applicability 
are for moderately collisional plasmas, as described by the Landau- Fokker- Planck equa- 
tion, the method is expected to be applicable also to collision processes described by the 
Boltzmann equation. This scheme is similar to the previously discussed velocity-based 
scheme [R. Caflisch, et. al, Multiscale Modeling & Simulation 7, 865, (2008)] and the 
scattering-angle-based scheme [A.M. Dimits, et. al, Bull. APS 55, no. 15 (2010, Abstract: 
XP9. 00006)], but with a firmer theoretical basis and without the inherent limitation to 
the Landau- Fokker-Planck case. It gives a significant performance improvement (e.g., 
error for a given computational effort) over the velocity-based scheme. These features 
are achieved by assigning passive scalars to each simulated particle and tracking their 
evolution through collisions. The method permits a detailed error analysis that is con- 
firmed by numerical results. The tests performed are for the evolution from anisotropic 
Maxwellian and a bump-on-tail distribution. 
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1. Introduction 

Any model equation of a system of kinetic particles necessarily contains a degree of 
error. The challenge, in the applied mathematics sense, is to control this error. For 
example, the Vlasov equation's errors are bounded to the degree that collisions may be 
ignored; the Euler, Navier-Stokes and Braginskii equations, to the degree that collisions 
dominate; the kinetic MHD equations, to the degree the gyrofrequency exceeds other 
frequencies. Each of these may be regarded as a perturbation expansion of the Boltzmann 
or Landau-Fokker-Planck (henceforth LFP) equation. 

Each of these expansions has the benefit of considerably simplifying the numerical 
simulation of the system in question. Of obvious interest is the extension of such compu- 
tational simplifications to regimes in which these expansions are not directly applicable. 
In one such regime - that in which collisions are important to the dynamics but not so 
dominant as to permit a fluid closure - the idea of combining Monte Carlo methods with 
fluid solvers to get accurate and efficient simulations has gained popularity in the last 



15 years or so [6J, [12|, |15|, ll7|, [18|, |20j, |24], |28|, |32|, with applications to both plasmas and 



rarefied gases. While useful in practice, the modeling considerations inherent in many of 
these schemes prevent a mathematical account of the size and scaling of the associated 
errors. 
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In this paper, we take a step toward a scheme applicable to this regime whose errors 
may be bounded in the same sense as in the perturbation expansions above. The scheme 
we present is a modification of those developed in that is more mathematically 

justified than either and more accurate than the scheme in [(J. While we treat the case 
of Coulomb collisions in a plasma - that is, the method generates approximate solutions 
to the LFP equation - an advantage of the present scheme over those in jgl, [ll| is that 
the core ideas may be applied to any elastic collision process described by the Boltzmann 
equation. In the present work, we treat the spatially homogeneous case, but the method 
is also not limited to this scenario. 

Moderately collisional plasmas appear in a variety of applications, including the toka- 



mak edge plasma |2ll. |29|. inertial confinement fusion [8], and counter-streaming astro- 



physical plasmas |5| . Generically, any system characterized by large variations in temper- 
ature and/or density is likely to feature a region in space where collisionality is moderate. 
Moreover, even in largely collisionless systems, collisional simulation may be necessary to 
model turbulence due to the small-scale structure developed jlj]. 

Typically, simulation of the full LFP equation is required when plasma collisionality 
is moderate. Particle-in-Cell (PIC) methods, in which the kinetic distribution / is repre- 
sented as an average over many simulated discrete particles, are often used jij]. Because 
of the complicated nature of the collision operator in the LFP equation, it is frequently 
replaced with a simplified model operator (e.g. [l|, Q, H] ) • Takizuka and Abe 3^1 (hence- 
forth TA) and Nanbu 27| have introduced binary collision algorithms that accelerate the 
simulation and have the additional advantage of approximating the LFP collision operator 
to O(At) |H, @|- There have also been advances in the Langevin equation-based treat- 
ment of the LFP collision operator (e.g. 25, 26[). Even so, the simulation of Coulomb 



collisions frequently represents a computational bottleneck in LFP simulations. Not only 
does the smallness of the collisional time scale restrict the time step size, but there may 
be multiple disparate collisional time scales in a single system jof. The necessity of cap- 
turing the shortest scale makes the observation of long time scale effects - sometimes the 
more important ones - very inefficient. We give an example that illustrates this point in 
section 2.2. 

The scheme we present here accelerates LFP collisional simulation by assigning passive 
scalars to each simulated particle which are evolved throughout the simulation. The 



scheme has commonalities with those presented in |13|, |15j, |24], [33| , but the development 
here is independent and less heuristic. The mathematical nature of the derivation makes 
it possible to perform a detailed error analysis of the scheme. 

The remainder of the paper is structured as follows: Section 2 presents background 
on the LFP equation and its asymptotic limits. Section 3 summarizes previous results 
on hybrid schemes of the type introduced in Section 4 motivates and outlines the 
steps in our new scheme. Section 5 details our methodology for tracking the values of 
the passive scalars assigned to each particle. Section 6 summarizes and presents an error 
analysis of the complete algorithm. Section 7 presents numerical results for two test initial 
conditions: a slightly anisotropic Maxwellian and a bump-on-tail distribution. Finally, 
Section 8 presents conclusions and indicates directions for future work. Some details are 
left to appendices. 
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2. Background - Kinetic Description of Plasma Systems 



Kinetic equations describing systems of many interacting particles take the form 

dtf + vd K f + a-d v f = C(fJ), (1) 

where /(x, v, t) is the particle number-density in position-velocity phase space at the 
position space coordinate x, velocity space coordinate v, and at time t, while a is a mean 
acceleration term that - in general - depends on /, and C(f,f) is a bilinear operator, 
depending only on the velocity variables, that accounts for inter-particle collisions. The 
method we propose is in principle applicable to any equation of this form. For simplicity, 
we restrict our attention to a single particle species, but the method is easily generalizable 
to multi-species systems. In this paper, we treat the special case of the LFP equation and 
indicate in the conclusion how the method may be applied to other collision operators. 



The LFP equation, originally derived by Landau 23], is obtained by setting 



8iielm 2 dv \ dv dvdv dv 



where H and G are known as the Rosenbluth potentials |3jj, which solve 

A V H = -Aug, A V G = 2H, (3) 

and e is the charge of an individual particle, m its mass, Eq the permittivity of free 
space, and log A is the well known Coulomb logarithm, with A - often called the plasma 
parameter - given by 

A = 4irn\ 3 D . (4) 

Here n is the particle number density in position space and A/j the Debye length, given 
by 

As = (5) 
V ne z 

where T is the system temperature. In a typical plasma, A 3> 1. 

The LFP equation is a standard kinetic model of plasma behavior when combined 
with the constitutive relations for a; namely, the Lorentz force and Maxwell's equations. 
Note that the collision operator fl2]) has an associated time scale 

(6^77/ \ 
A 2 2 3 lQ g A » ( 6 ) 

where v t is the magnitude of the typical relative velocity between two colliding particles. 



Analytic solutions of the LFP equation are known in only the simplest scenarios [30 
so one is forced to seek numerical solutions for situations of physical relevance. Even this 
effort is challenging due to the high dimensionality, non-linearity, and non-locality of the 
equation. This paper presents an accelerated Monte Carlo method for the treatment of 
the LFP collision operator (j2J). In describing the method, it is useful to first discuss two 
asymptotic regimes in which the problem is considerably simplified. 
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2.1. Asymptotic Limits 

We first consider the case in which tpp is asymptotically small compared to other 
time scales. Then, the system is highly collisional and permits a perturbation expansion 
based on the assumption that C(f,f) dominates ([T]). In this case, / is a Maxwellian 
distribution to leading order, given by 
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n m v — u 

/m(x ' v ' i) = (^F exp i — M^J- (7) 

where n, u, and T are functions of (x, t). In what follows, we denote a Maxwellian 
distribution by / m (v; n, u, T). Note also that for a Maxwellian, we may take v t = yj2T/m. 

Integrating ([T]) and using this expansion, one may obtain closed systems of PDEs for 
finitely many moments of /, the first few of which are defined by [?J 



n(x,t) = / fdv={f), (8) 

nu(x,i)= f vfdv=(vf), (9) 

If 1 

nT(x,t) = -m |v - u| 2 /dv = -m(|v - u| 2 /), (10) 

where we've also introduced (•) as a shorthand for velocity integration. For neutral fluids, 
the leading order term gives the Euler equations of fluid dynamics. Analogous equations 
are obtained for plasmas. The resulting equations are more amenable to numerical solu- 
tion than the LFP equation. 

In many plasma systems though, the opposite scaling holds: tpp is asymptotically 
large. The approximation C(f,f) = is made in (JTJ to obtain the Vlasov equation. 
Such systems are frequently called collisionless. PIC schemes are commonly used for 
solution of the Vlasov equation [i[ . 



2.2. Moderately Collisional Regime 

In the moderately collisional regime we study here, neither of the previous limits is 
valid. In such cases, PIC methods are frequently used with an additional step in which 
the simulated particles in a given spatial cell undergo collisions with each other. However, 
as already mentioned, the simulation of collisions frequently represents a computational 
bottleneck. Not only does tpp restrict the time step, but the collision time may vary 
within the system. 

Consider, for example, a so-called 'bump-on-tail' distribution given by 

/ = / TO (v;ni,0,7i) + / m (v; n 2 , u B , T 2 ) (11) 

with ub S> maxj y^ZT- The time scale tpp for intra-Maxwellian collisions is much shorter 
than that for inter-Maxwellian collisions, so the former dominate the computational effort 
but don't change the distribution. This makes direct Monte Carlo methods very inefficient 
for this and similar problems. 
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3. Hybrid Fluid-Monte Carlo Schemes 

Hybrid schemes of the type we consider arise from a splitting of the distribution / 
into / = ju + fk-, where fu is some initially Maxwellian distribution satisfying f^ < f. 
If f M and f k satisfy 

d t f M = C(f M , f M ) + C(f M , f k ) + S (12) 
d t fk = C(f k ,f k ) + C(f k J M )-S (13) 

where S is some arbitrary function of (v,t), then the spatially homogeneous version of 
flU is satisfied by / = f M + f k . 

Since we choose our splitting such that /m is initially Maxwellian, C(/m,/m) = 
initially. If /m happens to remain close to a Maxwellian throughout the evolution, we 
are justified in ignoring C(/m,/m) completely. In this way, this splitting generates a 
computational savings over PIC by avoiding the necessity of simulating collisions between 
particles within fu- If I'm constitutes a large fraction of the system's mass, then this 
saving will be significant. The scheme also gains accuracy over a pure fluid scheme, 
because such schemes treat only perturbative deviations from a Maxwellian, while the 
split kinetic system (I12|) -(13) does not assume f k is asymptotically small. 

The problem, then, becomes choosing S in such a way that fu remains - at least 
approximately - a Maxwellian for all future times, while at the same time maximizing 
the fraction of the system's mass residing in fu- The closer S keeps /m to a Maxwellian, 
the more accurate the scheme. The more positive S is, the faster the scheme. 

Previous efforts resort to model equations for S 0, EH . These models are more easily 
summarized if we rewrite 

S = f k r T ~ fur D (14) 

for some positive and representing distribution normalized v-dependent exchange 
rates into and out of /m, respectively. 

In the presence only of collisions, the system tends toward a Maxwellian distribution. 
The quantity tt represents the transfer of particles from f k to $m to reflect this. Increasing 
7Y makes the scheme more efficient but less accurate. The movement of kinetic particles 
into f\j will be referred to as thermalization. 

Similarly, collisions between particles in /m and f k drive fM away from its current 
equilibrium state. The quantity m represents the transfer of particles from f^ to f k 
to reflect this. Decreasing m makes the scheme more efficient but less accurate. The 
movement of particles from }m into f k will be referred to as dethermalization. 

We shall refer to any method that prescribes and To - and therefore the number 
and variety of particles to be thermalized and dethermalized at each time step - as a 
thermalization scheme. The derivation and testing of an improved thermalization scheme 
is the subject of this paper. 

We first discuss previous thermalization schemes used in hybrid methods of the type 
described here. Since these methods require finite time steps, we choose to write our 
statements in terms of 

p T = r T At, p D = r D At, (15) 

the probabilities of a given particle being thermalized or dethermalized in a given time 
step of length At. The variables on which px and po depend characterize previous 
thermalization schemes. 
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3.1. Velocity-based Schemes 

Introduced in ((jj], velocity-based schemes have px, Pd dependent only on \w M — v p |, 
where v p is the particle's present velocity and Um is the mean velocity of /m- Pt is a 
decreasing function of this quantity, while po is increasing. 

This scheme is intuitively sensible, but has numerous drawbacks. Firstly, there are 
many choices for p? and pjj, and it is unclear if an optimal choice even exists. Just as 
serious, we claim, is the conflation of "similar velocity" with "many collisions" . It is true 
that if a particle undergoes many collisions with particles from a given Maxwellian, its 
mean velocity will tend toward that of the Maxwellian. However, the converse is most 
certainly false. To illustrate this, consider the following initial distribution 

/ = /m(v; ni, 0, T) + / m (v; n 2 , 0, eT) (16) 

for e < 1 and n%<n\. A hybrid method might divide this distribution into a Maxwellian 
component (/m) given by the first Maxwellian and kinetic component (/&) given by the 
second Maxwellian. 

If a velocity based scheme is to be efficient for the example (fl6|) . it should immediately 
thermalize every kinetic particle, because their velocities are very near the center of Jm- 
On the other hand, this is not possible since if a velocity based scheme is to be accurate, 
Pt must be relatively small, even at /m's mean velocity, thereby sacrificing efficiency for 
other initial conditions. We conclude that the velocity of a particle alone is not enough 
information to decide whether or not it should be thermalized. 



3.2. Scattering Angle-based Schemes 

The thermalization scheme developed by Dimits et. al JJJ is such that pt depends 
only on 9, the scattering angle the particle subtended in its most recent collision, and one 
additional (overall multiplier) parameter. In its current form, this scheme sets po = 0. 
Pt is typically an increasing function of 9. 

This scheme is intuitively sensible for the case of Coulomb collisions in which small an- 
gle collisions dominate the dynamics. The applicability of this scheme to other potentials, 
for which small angle collisions do not dominate, is questionable. 

Like velocity based schemes, scattering angle schemes use only velocity information 
from the most recent time step in making decisions about thermalization. We argue that 
it is desirable for a scheme to make use of additional variables, which better capture the 
long-term collisional history of the kinetic particles. 



4. Paradigm and Theoretical Background 

As discussed in the previous section, we claim that a particle's velocity is not enough 
information to determine whether it should be thermalized, nor is any information depen- 
dent on only the most recent time step. We claim that one should look at the distribution 
of velocities that particle might have had, given its collisional history. This point merits 
elaboration. 

A Monte Carlo scheme represents / as a sum of particles with known velocities. Each 
particle undergoes a sequence of random collisions throughout the simulation, so that 
after any given number of time steps, the velocity of a single simulated particle may be 
regarded as a random variable. Let us denote by fj(v, t) the probability density function 
of the velocity v of the j'th simulation particle's velocity at time t, rescaled so that the 
total mass of the f/s matches that of /. Notice that this is initially a delta function at 
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the particle's designated velocity. In a Monte Carlo scheme, / is realized - conceptually 
- as 



where the sum is taken over all simulated particles. 
The analogous equation for the hybrid scheme is 



f = fu + J2fr 



(17) 



Moreover, we may apply a splitting analogous to that in ( 1T2"]) -(13): 

dtfM = C(f M J M ) + Y^C(f M , 

j 
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%fi = 'E l C(f j ,f i ) + C(f J ,f M )-S j , 



(19) 
(20) 



where the indices j and % run over all the simulated particles. We will often write /& = 

Z)j fj- 

In this framework, once (I19|) -(2Q) are discretized with time step At, the thermalization 
of the jth particle amounts to setting 



f M {t + At)=U M {f M {t) + f, 



3)1 



(21) 



where Il^f is the projection operator onto a Maxwellian - that is, IIj^J is a Maxwellian 
with the same (n, u, T) as /. The goal of the scheme we propose is to thermalize particles 
in such a way as to introduce as little error as possible into the overall scheme. That is, 
the decision to thermalize a particle should have little effect on the overall distribution. 
This is achieved by thermalizing the jth particle only if 

1 



\\f M + fj-UM(fM + f; 



3 JUL} 



< £ 



(22) 



for some e > 0, where rij = (fj)- The choice of norm is somewhat arbitrary, but L 1 will 
prove convenient later. By the triangle inequality, (122]) is implied by 

fj 









/at _ fj 


+ 




LI 


V n M J 



' M 



< TljE, 



(23) 



where ju = ( n j/nM)fM- The second norm is small if fj ~ /m, which it must be for the 
first norm to be small. It's even smaller if rij um, which is the case in the parameter 
regimes we consider. We therefore find it sufficient to enforce 

1 : u 



f 



M 



< e 



(24) 



as a condition for thermalization. 

We propose to thermalize the jth particle whenever §M§ is satisfied and to dether- 
malize it when (124]) is violated. In order to implement this, we need a way of computing 
the L 1 norm in (124]) . This is not a simple task. We instead compute a quantity called 
relative entropy, which bounds the L 1 norm from above, and thermalize particles when 
this quantity is sufficiently small and dethermalize them when it is sufficiently large. Rel- 
ative entropy has the added advantage of evolving monotonically through the action of 
collisions with a Maxwellian background, as we show below. 
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4-1. Relative Entropy 

Define the relative entropy (sometimes called the Kullback-Leibler divergence) %{f\g) 

by 

H(f\g)= [ ]ag(£)fdv. (25) 

for any two non-negative functions /, g depending on v. This quantity has its origins 
in information theory, where it is used in the context of coding information sources jij. 
Here, we need only cite four properties: for any non-negative / and g satisfying (/) = (g), 



U{f\g) > (26) 

H(f\g) = iff f = g (27) 

\\f-9\\l<2(f) 2 n(f\g) (28) 

C{f,f m )\og(V) dv<0 (29) 



where f m is a Maxwellian, f m = ((/)/ (f m ))fm,J = // (/) and similarly for g. 

The first two properties are standard (e.g. [9j). The third is a straightforward gener 



alization of the CKP inequality jlOL |22| to distributions with non-unit mass. The fourth 
is a modification of Boltzmann's if-theorem whose proof we present in appendix A. 
The inequality (125]) allows us to impose fpH|) by bounding "H(/j |/m), since fl28|) implies 

U{f 3 \] M ) < £ - =► - / M - £ „ < e. (30) 
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The result ( 129]) states that the role of the term C(fj, Jm) i n (19) is to drive fj irreversibly 
toward fu- Moreover, (129]) shows that collisions between / and fu tighten the L 1 
bound monotonically in time. We have not shown that C(fj,fnj) drives fj toward /m 
at any particular rate. However, numerical experiments in section 7 suggest the rate is 
comparable to tp P . The final thermalization criterion is 

'Hi f, f \i) < H r . (31) 

where "H c is a parameter of the scheme. We specify the scale of this parameter in section 
6.1. 

4-2. Idea for Thermalization Scheme 

This leads us to a more precise formulation of the thermalization scheme previously 
proposed: To the jth kinetic particle fj, assign a passive scalar representing its relative 
entropy, H(/j|/m)- Whenever it undergoes a collision in a given time-step with a particle 
whose distribution is given by fa, we evolve this passive scalar in a way that is consistent 
with (20). If at the end of any time step the particle's relative entropy has dipped below 
the threshold T-L c , we thermalize it. 

Similarly, whenever a particle's velocity must be sampled from the Maxwellian compo- 
nent in order to perform a collision, we assign it a relative entropy of zero (corresponding 
to fj = /m), and evolve it one time step according to the same kinetic equation. If at 
the end of the time step, the particle's relative entropy exceeds T-L c , we dethermalize it, 
and it carries with it this relative entropy value. 

It remains to specify the details of evolving the relative entropy of fj according to 
collisional terms in (20). The following section is devoted to doing this for the case of 
Coulomb interactions. 
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5. Approximating Relative Entropy 

Toward evaluating "H(/ j |/m) ) we n °t e from (19) that its rate of change due to collisions 
with f M during a single time step is 

{W)m = [ log (4-) C (fv fu) dv, (32) 
J® 3 V/m/ 

where the subscript M indicates that we're only treating collisions with Jm- In 5.3, we 
show how to extend this treatment to the other terms in (20). 

The right hand side of fl32l is, in general, difficult to evaluate. In particular, evaluating 
the right side exactly requires knowledge of the full distribution fj for each kinetic particle, 
which is not computationally feasible. We approximate the integral by approximating fj 



by a finite-moment truncation of a tensor expansion [19|. The more moments we keep, 
the better the approximation of fj, but the more quantities we have to evolve at each 
time step for each kinetic particle. We make the compromise of keeping the standard five 
moments defined in (jHJ)- (10), corresponding to the assumption that fj is a Maxwellian. 
This is justified both early in the simulation - when fj m 5 3 (y — Uj) - and late - when 

fj ~ /m- 

We will say that /m has temperature Tm and mean velocity um, while fj has mean 
velocity Uj and temperature Tj. Algebraic manipulation of ( 12 5 p reveals that 

nfj\f M ) = I ( 7 I^ + lo g (^-)) + (33) 



2 \ T M \Tj J J 2T M 

where m is the common mass of all the particles under consideration and UjM = Uj — u« . 
To specify the relative entropy in this case, it is thus enough to specify Uj and Tj, so 
instead of having to compute the integral in ([3"2"j) . we can work with the comparatively 
simple collisional rates of change of mean velocity and temperature: 

(d t Uj) M = vC F p(fj, f M ) dv, (34) 

{d t Tj) M = \m |v - Uj\ 2 C FP {fj, f M ) dv. (35) 

This approach has the additional advantage of avoiding the direct approximation of d{K, 
which can be arbitrarily large, as indeed can "H itself (consider Tj — > in (I33p ). Deriva- 
tives of Uj and Tj are much more well behaved. 

Some results on the integrals in (1341) and fl35|) are derived in [lij], and some others are 



attributed to Decoster in [20|, [24J]. Here, we make use of both results as well as deriving 
some that are new to the best knowledge of the authors. Because the derivations are 
lengthy and the results partially known, we leave the details to appendix B and present 
only the results here. 

Define FjM as the right hand side of (1341) . and 2W}m/3 the right side of (1351) . When 
Tj Tm, we show in appendix B that 



7 ~F<5 4 7 n M UjM 

in u tM u jM 



jM dx en \ u jM) 



(36) 



W jM « W} M = 2 -^L^l. (37) 

J mv tM Uj M 
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where U,m = ^ju/vtM, with v t M = a/ 2Tm / m denoting the thermal velocity of fu- For 
the definition of the constant 7, see (IB.3j) . Equation ( 136]) is equivalent to the analogous 



result in [20|,|24] 



When UjM with v t j denoting the thermal velocity of fj, we show in appendix 

1 that 

(38) 
(39) 



jM 



w. 



jM 



w m 



~ V jM 


1 

TjM 









(T M - Tj 



ma 



•jM 



where 



TjM 



3v^Fm 2 (v§ + v 2 tM ) 



3/2 
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7 n Af 



Note that f[3"81) is equivalent to the analogous result in [141 ]. while 
of the analogous result in same. 



(40) 

is a generalization 



5.1. Toward a uniformly valid approximation of relative entropy 

While the above asymptotic expressions are interesting in their own right, we require 
expressions that are uniformly valid throughout parameter space. The integrals in (13~4"|) 
and (I35p can be calculated numerically, but the inline evaluation of mult i- dimensional in- 
tegrals is computationally prohibitive in this context. It is possible to reduce the problem 
to one dimensional integrals dependent on only two non-dimensional parameters (see ap- 
pendix C), from which a look-up table can be generated. However, numerical experiments 
presented in section 7.1 show that the majority of the error in our relative entropy estima- 
tion comes from the assumption that fj is Maxwellian, rather than from the asymptotic 
assumptions used to derive fl36|) - fl39|) . We will therefore find it satisfactory to set 



(d t u 



j)M 



M 



2 
3 



Tp<5 

b jM ■ 
Tpm 

b jM ■ 

w 3 5 M 
w. 



UjM > CtV tj 

UjM < av tj 



Uj M > OtV t j 

jM ■ u jM < no, 



(41) 



(42) 



where a G (0,1). This clearly recovers the correct behavior when UjM *C ihj. To see 
the recovery of the other limit, we rely on numerical experiments to confirm that when 
UjM 3> v t j, it is also the case that Tj T M , making fl3"6"|) - fl3"T|) valid. All tests indicate 
that the properties of the scheme are not sensitive to the choice of a. In all results that 
follow, we use a = 0.9. 



5.2. Monte Carlo Implementation 

Thus far, we have discussed the role of the term C(fj, /m) in (fl9|) -(20) in the evolution 
of Tj and u,-. There are two other collisional effects that must be taken into account. 

Firstly, we must also evolve Tj and u, through collisions with other kinetic particles. 
That is, we must account for the term C(fj, fi). This is done by retaining our assump- 
tion that each of the particle distributions is a Maxwellian. Then, the difference between 
treating collisions with /m and with fi is merely a matter of changing parameters, since 
both satisfy the assumptions in previous subsections. For instance, if we were to rewrite 
to give the rate of change of Tj due to collisions with the ith particle, we would have 



W" 1 = — 

iji 



Tj 



mu 



.1 1 



(43) 
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where Uji = Uj — Uj, and 

a^ M^lg . (44) 

16 7?2j 

The expression (I43|) is then valid when Uji <C v t j. Analogous changes are made to (I36|) - 
fl38|) . and when the jth kinetic particle collides with the ith kinetic particle, we evolve Uj 
and Tj according to 

<*»'>< = \ $ : :z < r < 45 > 

(#Z& = MS : UjM - ^ (46) 
v 3,1 3 \W% : UjM < av tj v ; 

where the % subscript denotes change due to ith particle. 

Secondly, we have so far only looked at the rate of change in Uj and Tj for one of the 
two collision partners. The collision also affects the moments of the other distribution. 
These effects are captured by recalling that collisions conserve momentum and energy. 
That is, 

(v [C(f,g) + C{g, /)]} = 0, (47) 
(v 2 [C(f,g) + C(g,f)}) = 0, (48) 

for any / and g. Therefore, we have 

Fjt = (vC FP {fi, fi)) = -(vC FP (U fj)) = -Ftf> ( 49 ) 

and 

2m~ l W ji = 2m- 1 W ij - 2u j{ ■ ¥ {j . (50) 

In this way, our treatment of the terms C(fj, J'm) and C(fj,fi) implicitly generates the 
analogous formulas for C(/m, fj) and C(fi, fj). 



6. Algorithm Summary and Error Scalings 

Denote the number of simulated particles constituting by Nk, and that constituting 
fu by N m . To generate a numerical solution of the spatially homogeneous LFP equation, 
at each time step the following substeps are taken. 

1. Determine the number of collisions of each type to perform. In the algorithms 
of Nanbu and TA, each particle undergoes exactly one collision with a randomly 
selected partner: 

(a) Nkk = N k 2 /2(Nk+N m ) is the number of collisions between two kinetic particles. 

(b) N m k = N m Nk/ (N m + Nk) is the number of collisions between a kinetic and a 
Maxwellian particle. 

(c) Collisions between Maxwellian particles do not affect the distribution and are 
thus not simulated. 

2. Perform collisions: 

(a) Randomly select Nkk kinetic particles, and assign each a partner from among 
those not already selected. Perform collisions between each pair by altering 
their actual velocities according to the collision algorithm as well as their mean 
velocities and temperatures according to flj5j) and (JIB")) , respectively. 
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(b) For each kinetic particle unused in (a), generate a particle with velocity sam- 
pled from /m, and with mean velocity and temperature equal to those of Jm- 
Perform a collision between this and the kinetic particle as in (a), but evolve 
the mean velocity and temperature of each particle according to (|4"Tj) and (l4"2l . 
respectively 

3 . T hermalizat ion /Det hermalizat ion : 

(a) For each kinetic particle, compute from Uj, Tj using ( |33i) . If this 
quantity is less than % c , remove the particle from the kinetic component, 
increment N m and decrement iV fc . 

(b) For each particle sampled from f M , compute "H(/j |/m) as in (a). If this number 
exceeds H c , add the particle to the kinetic component, decrement N m and 
increment N^. 

4. Enforce conservation: 

(a) Adjust um so that the total momentum in the system is the same as before 
the collisions. 

(b) Adjust Tm so that the total energy in the system is the same as before the 
collisions. 

If there exists spatial dependence in the problem, fluid equations which adjust tim, Um, 
and T M must also be evolved. 

6.1. Choice of% c 

The algorithm outlined above has one free parameter: % c , the value of relative entropy 
below which a particle is thermalized and above which it is dethermalized. We now present 
an argument that specifies the scale of this quantity, although not its precise value. 

Notice that, unlike the thermalization process, the dethermalization process occurs 
over a single time step. That is, if the scheme doesn't dethermalize a given particle, 
that particle is inserted back into the Maxwellian and its interaction with the kinetic 
component of the scheme is forgotten. Accuracy demands that this does not happen to 
every particle sampled from the Maxwellian, so that some particles get dethermalized. 
Thus, T-L c should not be much larger than the change in relative entropy experienced by 
a Maxwellian particle in a single time step. 

Similarly, efficiency demands that not every sampled particle be dethermalized. There- 
fore, He should not be much smaller than the aforementioned change in relative entropy. 
We conclude that T-L c should be comparable to the typical change in relative entropy 
experienced by a particle sampled from the Maxwellian. It remains only to specify this 
scale. 

Denote the change in temperature a Maxwellian particle undergoes during a collision 
with a kinetic particle by AT and the change in u by Am. If we assume AT ' /Tm <C 1, 
then to leading order in ( 133|) . the change in relative entropy is 




(51) 



From ( |39|) and (|50|) . we can estimate AT. 



AT 



8 



7 



T M At 



(52) 
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Similarly, by using ( 138]) . we can estimate 



3 7 

Am « - _ ' uAt (53) 

By rewriting the above in terms of t FP using flBJ) and f)B.3j) . then plugging into f lBTj) . 
we have 

12 / 3 / u s 




H.„AH» T ^l + 5S j^— Ml-) »3^|-| . (54) 

In the last approximation, we've assumed that u is not so large as to make the second 
term significant, which is valid so long asu< 9v t M- In general, we set 

At x 



H c = c^-j (55) 

for some c we choose. Unless otherwise specified, we use c = 12.9 henceforth. This choice 
makes the errors more visible, although the fidelity of the results is not sensitive to this 
choice. 

This scaling for T-L c implies that the (de-)thermalization error scales like At, since by 
f[3"0j) we have 

m.1) fit) < O =► 1 \f 3 - hi < O . (56) 



In numerical experiments to follow, we find that the accuracy and efficiency of the 
scheme are each insensitive to changes in T-L c within a factor of 10 to 100 around the value 
in (GET 



6.2. Error Analysis 

For a spatially homogeneous PIC scheme using the collision algorithm of TA or Nanbu, 
we have 

Wftme - fpic\\ Ll = 0(At/t FP ) + 0(N~^ 2 ), (57) 

where f trU e is the analytic solution and fpic its approximation by PIC with N particles 
and time step At. 

A hybrid method which never thermalizes particles and dethermalizes every particle 
sampled from the Maxwellian component reproduces the result of the PIC scheme in 
expectation. Moreover, because each simulated particle carries equal weight, we have 
rij = n/N. Thus, each thermalization or failed dethermalization event introduces an 
0(At/tppN) error by f)56p . There are O(N) failed dethermalizations at each time step, 
and 0(1 /At) total time steps, so there are 0(N/At) failed dethermalizations over the 
course of the simulation. 

Contrast this with the scaling of the number of thermalization events. This number 
may be said to be J z ('H c )0(N), where J 7 is the fraction of the simulated particles that 
are thermalized during a simulation, which is clearly an increasing function of T-L c . Since, 
in addition, T-L c is an increasing function of At, we may say that J 7 = 0(At") for some 
j3 > 0. This gives the overall scaling 



11/ 



PIC 



Mli = o (^) (o g) + o (n*t>) = o (#) + o (^) . 

(58) 
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Combining (15 7j) and fl58|) gives the error scaling for the hybrid scheme proposed here: 



/ At\ , /Af 1+/3 \ 

Wftrue ~ fhy bnd \\ Ll = O f — J + O {N^' 2 ) + O (tpp) + O I -j— J 



(59) 



where the first term is the finite time step error, the second the sampling error, the third 
the (failed) dethermalization error, and the fourth the thermalization error. 

The derivation of (1391 has the following weakness: it assumes that we know the actual 
value of the relative entropy for each fj, when in fact we only have an estimate of it based 
on the assumption that fj is Maxwellian. The derivation of (I59p remains unchanged if 
we have 

n true < m est (qq) 

for some k > 0, where the superscript true indicates the actual relative entropy for a 
given particle, and est indicates our estimate of that quantity 

However, for any 7-L est computed using a finite moment truncation of fj, there exist 
pathological cases in which l-i est = while % true is strictly positive, so ( 16 Op cannot 
hold in general. However, for many problems these pathological conditions may not be 
achieved, giving hope that (l59|) may be realized. Indeed, results in the following section 
are consistent with the scaling relations presented here. 

6.3. Dethermalization Error 

Since failed dethermalizations are shown in (|59|) to be the dominant error source, it 
warrants more discussion. In particular, having seen that this error doesn't scale with 
At or N, we seek to understand what sets the size of this error. 

By (1251) and (1541) . we have a bound on the error incurred by each failed dethermaliza- 
tion event, which we denote by e e : 

i—n At 

^ £ ^nT7 p (61) 

for some constant c. The growth rate of this error, denoted by e t , is thus given by 

e t = e e TZ, (62) 

where 1Z is the number of failed dethermalization events per unit time, which may be 
written as 

K = — — 63 

n At y ' 

when rik/n is small. Combining ( 16T]) - (I63|) . we have 

e t \f2cnu , JN 

— < -, 64 

n tpp n 

where we divide through by n so that the right side may be thought of as a fractional 
error. 

The tpp above is the characteristic time for collisions between and fk, which will 
change throughout the simulation. We find it instructive to express the same statement 
in terms of the characteristic time for collisions within /m, denoted by tp P . 



£t V2c n k ( v tM x 
n ~ tf p n \u kM 



< Tm~— ( — ) » ( 65 ) 
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where u^m is the characteristic relative velocity between a particle sampled from and 
one from ju- 

The inequality fl65|) highlights two ways in which the dethermalization error can be 
made small: the kinetic component can be small - i.e. <C n - and/or the kinetic 
component can have velocity very different from the maxwellian component - u^m ^> VtM- 
Moreover, the case in which neither of these conditions is realized is short lived, as a 
kinetic component with UkM ^ VtM will be rapidly thermalized. 

We compare results from numerical tests to the predictions of ( |65|) in the following 
section. 

7. Numerical Results 

We perform two types of numerical test. First, we check that the approximations ( 1331) . 
( |4T|) . and ( |42|) capture the actual evolution of relative entropy through collisions with a 
Maxwellian background. Second, we test the entire algorithm's accuracy and efficiency 
against pure PIC simulations in two test cases: the relaxation of a slightly anisotropic 
Maxwellian and a bump-on-tail distribution. 

7.1. Testing Relative Entropy Approximations 

We are interested in the following problem: given a distribution /(v,t) solving 

d t f = C FP (fJm), /(v,0) = 5(v-v o ), (66) 

where f m is Maxwellian and constant in time, what is the time evolution of rL{f\f m )l 

We attack the problem in three different ways. First, we represent / by a single test 
particle with initial velocity v , and evolve its velocity according to the algorithm of TA, 
where each collision partner has velocity sampled from f m . The test particle's velocity is 
then a random variable v p (t) whose distribution is given by f(y,t). By simulating the 
evolution of v p repeatedly, we may generate - at each t - a histogram that, by definition, 
approximates / at that time. We then evaluate rl(f\f m ) by direct numerical integration 
of (1251) . This will be thought of as the true value of the relative entropy, as it is subject 
only to the errors in the Monte Carlo scheme and the numerical integration, each of which 
can be made arbitrarily small. Results from this method are plotted in solid blue in fig. 
1. 

Second, we evolve the mean velocity and temperature of / according to the collisional 
moments F and W under the assumption that / is a Maxwellian. F and W are numeri- 
cally integrated (refer to (IC.5P and (IC.lip in appendix C) to find the rates of change of 
u and T, which are then evolved by forward Euler and plugged into fl33l) to find ri. This 
method is designed to test the validity of the assumption that / is Maxwellian, since this 
is the only assumption made here that was not used in the direct simulation of collisions 
outlined above. Results from this method are plotted in dash-dot black in fig. 1. 

Third, we evolve u and T according to the numerical solution of the ODEs (|4ip and 
fl4"2j) . then evaluate ri using (1531) . This tests the validity of our asymptotic expressions, 
since this is the only assumption used here but not in the previous method. Results from 
this method are plotted in dashed red in fig. 1. 

Fig. 1 shows the values of u, T, and "H^I/m) computed from each of the three meth- 
ods above for three different values of Vq. From this, we conclude that the approximations 
(1411) . (I42p are a satisfactory framework for the approximation of relative entropy 
evolution. In particular, the approximation captures the monotonicity and overall rate 
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Figure 1: Comparison of scaled temperature (T), mean velocity (u), and relative entropy (H) 
of the distribution of a test particle in a Maxwellian background computed by three different 
methods: direct Monte Carlo collision simulation (solid blue), evolution via numerical inte- 
gration of F and W (dash-dot black), and evolution according to the asymptotic expressions 
([IT]) and (fl2|) (dashed red). The left-most column has Vq = 0, the middle vq = vtM, and the 
right-most vq = 2vtM 

of relative entropy decay, and is especially accurate when the relative entropy is small, 
which is of particular importance for our application. 

Increased accuracy in the evolution of T and u may be obtained through numerical 
integration of f]B.5j) . (IB. 14[) . which may be of independent interest, but this does little 
to improve the relative entropy approximation, and in fact even degrades the quality 
of the approximation in some regimes (see bottom right plot in fig. 1). The numerical 
evaluation also greatly increases the complexity of the scheme in practice, so we use the 
asymptotic approximations (1331) . (14T|) . and (1421) in all numerical tests that follow. 

1.2. Full Numerical Tests 

We now present numerical results for the full algorithm outlined in section 6 when 
applied to two test problems: a slightly anisotropic Maxwellian, and a bump-on-tail 
distribution. 

As the goal of this algorithm is to accelerate the simulation of collisions, a discussion of 
computational cost is warranted. To leading order, the computational cost of collisional 
simulation is proportional to the number of collisions simulated during the scheme, which 
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is proportional to the number of simulated particles, averaged over the full run. Therefore, 
a hybrid scheme of the type we discuss is faster than a pure PIC scheme by a factor of 
roughly 

- /N\ / N m + N„ \ 

where the angle brackets now represent an average over all time steps. This is the measure 
of efficiency we use when testing schemes of the types outlined in sections 3.1 and 3.2. 

However, the entropy based scheme proposed here incurs an additional computational 
load for each collision due to the passive scalars that must be evolved. The cost of the 
simulation of any given collision is roughly proportional to the number of scalar quantities 
that must be evolved. This implies that the analogous efficiency measure for the entropy 
based scheme is 

- 3 / N m + N t \ 

S » = IT<\-^- /' (68) 

where d u is the number of components of u that we track, which will be problem de- 
pendent. In most cases, d u is the number of spatial dimensions in the problem since 
velocities in other dimensions are assumed to vanish on average. However, in problems 
with no spatial dependence but non-Maxwellian initial data, we will require d u ^ 0. In 
the bump-on-tail problem in 7.2.2, for instance, we set d u = 1. 

All of the following results were performed in a dimensionless formulation with m = 1, 
tp P = 5.348275, T = 0.05065776, and n = 0.1 (consistent with parameters in 

7.2.1. Two-temperature Maxwellian Relaxation 

We first test the fidelity of our implementation in a scenario with a known analytic 
solution. Consider the initial distribution 

nm 3/2 ( m(vl + Vy)\ ( mvl 



f ~ (2n)^TVTT5T eXP { IT 2(T + 5T) ) (69) 

with 8T <C T. In j^Ej], Trubnikov showed that the temperature difference 5T - to leading 
order - decays exponentially in time: 

5T(t) = 5T{0)e~ t/r , (70) 

with r given in Gaussian units by 



5 y/^T 3 / 2 

T 8v^Fne 4 logA' 1 ' 

In fig. 2, we compare this exact solution to the method of TA and the entropy- 
based hybrid method proposed here. We use <5T(0)/T = 1/10, At = tp P /20, and iV = 
1.024 x 10 6 with the other parameters as described above. We find agreement between 
all three solutions up to the level of statistical fluctuations in the numerical solutions. 

7.2.2. Accuracy and Efficiency Tests 

The bump-on-tail initial distribution we treat is given by 

f(t = 0) = / m (v; fa, 0, T) + / m (v; (1 - /3)n, u fc , T k ) (72) 
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Figure 2: The temperature anisotropy in a two temperature Maxwellian as a function of time. 
Computation compares the linearized analytic result (solid black), pure PIC using Takizuka-Abe 
(dash-dot red), and the hybrid scheme proposed herein (dashed blue). 

with (3 G (0, 1) corresponding to the fraction of the total mass in each Maxwellian. We 
set the initial Jm equal to the first Maxwellian term and equal to the second, and 
display plots for j3 = 0.9, u k = 2.83v t M^-, and T k = 10~ 4 . 

In fig. 3, we plot a time series of the hybrid solution compared to the Takizuka-Abe 
solution using At = tp P /20 and N = 256,000 total particles. The plots show excellent 
qualitative agreement between the PIC solution and the hybrid solution with Sh ~ 10. 

In fig. 4, we compare the efficiency and accuracy of the entropy-based scheme proposed 
here to the scattering angle-based scheme outlined in section 3.2. Analogous results for 
the velocity-based scheme described in section 3.1 may be found in [(J, showing inferior 
performance compared to both schemes tested here. For each scheme, the accuracy is 
measured by 

ptmax 

Tacc = ~ / \\fpw(t) ~ fhybrid(t)\\ L i dt, (73) 

"maxTl Jo 

the same measure used in 0. We set t max = lltpp to capture most of the progress 
toward equilibrium shown in fig. 3, although the results we present are not sensitive to 
this choice. 

For the scattering angle-based scheme, we set 

Pt = min |/csin -, 1 j> , (74) 

where 9 is the scattering angle in the two-particle center of mass frame, and vary k to 
change the efficiency - S - of the scheme. For the entropy based scheme, we vary the 
efficiency Sh by varying H c about the value prescribed in ( 15"4"|) . 
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Figure 3: Six time snapshots of the velocity distribution in the x dimension for the bump-on-tail 
problem. PIC solution: solid red; hybrid solution: black "x" . 

We test each scheme in two cases. In the first, no particle is ever dethermalized. In the 
second, we dethermalize particles according to po = Pt/2 for the scattering angle based 
scheme and as described in section 6 for the entropy based scheme. This is intended to 
test whether collisionally driven dethermalization plays a significant role in the evolution 
of the distribution and how efficiently each scheme handles the dethermalization process. 
Results for At = t^ P /20 and N = 256, 000 are shown in fig. 4. The sampling error is 
estimated at 0.015 by comparing multiple independent PIC simulations. 

We see immediately the improved accuracy effected by the entropy-based scheme for 
a fixed number of simulated collisions by comparing the green and red curves. However, 
the additional computational load incurred by the entropy scheme effectively cancels 
this gain in the absence of dethermalization, as seen by comparing the blue and red 
curves. However, in the presence of dethermalization the advantage is restored - seen by 
comparing magenta and black. 

Moreover, we notice that adding dethermalization has no effect on the efficiency of 
the entropy scheme (compare blue and magenta curves in fig. 4), while it degrades 
that of the scattering angle scheme (compare red and black curves). Dethermalization 
slows the scattering-angle scheme because dethermalized particles are slow to be re- 
thermalized. In the entropy scheme, on the other hand, recently dethermalized particles 
are re-thermalized quickly because they have particle temperature and mean velocity very 
close to those of the Maxwellian. We discuss the consequences of these observations in 
more detail in section 8. 

7.2.3. Convergence Study 

In fig. 5 we present a numerical convergence study comparing the various sources of 
error in (I59p . again using the bump-on-tail initial distribution. In an effort to isolate the 
systematic errors, we increase the number of simulated particles to N = 2.5398 x 10 6 . 
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Figure 4-' A comparison of accuracy T acc and efficiency S, Sh for four different realizations of a 
hybrid scheme: the entropy based and scattering angle-based schemes, each with and without 
dethermalization. The entropy scheme is plotted against both S (green '+') and Sh (blue and 
magenta dots). The scattering angle scheme is plotted both with thermalization (dashed black) 
and without (solid red). Comparing the green '+' curve to the solid red and dashed black 
compares the entropy and scattering angle schemes when the number of simulated particles is 
equal. Comparing the dotted blue to the solid red compares the two schemes when the total 
computational load is equal and when there is no dethermalization. Comparing the dotted 
magenta to the dashed black compares the two schemes for equal computational load when 
thermalization is present. A sixth possible curve - Entropy scheme without dethermalization 
vs. S - is not plotted, because it falls directly on top of the green '+' curve. 



We use time steps At = tp P 2~ k for k = 2,..., 7. The sampling error is of course not 
completely eliminated, and is thus subtracted from each curve in fig 5. For each curve, 
the sampling error is estimated by comparing multiple independent simulations. 

The red, unmarked curve in fig. 5 shows the errors between PIC schemes at the 
various time steps and the PIC scheme at the finest time. It is consistent with the 
expected O(At) convergence. The black, x- marked curve shows errors for the entropy- 
based hybrid scheme with no thermalization (we simply skip step 3a in section 6), i.e. it 
shows the dethermalization error. As expected, this curve is asymptotically constant, and 
begins to show the time-stepping error only to the right of the plot when time-stepping 
error becomes comparable to the dethermalization error. 

The blue, square-marked curve shows errors for the entropy based hybrid scheme 
exactly as summarized in section 6. The difference between the blue and black curves 
is the thermalization error, shown in the green, triangle-marked curve, and is found to 
scale like o(At), as predicted in (l59i) . All hybrid simulations presented in the plot have 
Sh > 5.9, with the full hybrid simulations (blue) having efficiency as high as Sh ~ 8.6. 
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Figure 5: A comparison of the errors incurred by PIC and hybrid schemes. 

We notice that the time stepping error appears to be smaller by a constant factor 
for the hybrid scheme as compared to the PIC simulations - i.e. at the right edge of 
fig. 5, the blue, square marked curve lies slightly below the red, unmarked curve. We 
hypothesize that this is because only a small portion of the distribution is subject to the 
time stepping error in the hybrid scheme, while the whole of the distribution is subject 
to it for PIC schemes. 

7.2.4- Dethermalization Error Study 

Lastly, we test some of the predictions of (IrJoT) using the bump-on-tail distribution. We 
investigate the rate of error generation by plotting the L l difference between the hybrid 
and PIC solutions as a function of time, using only dethermalization. We vary Uk in (1721) . 
which is analogous to uum in fl65l) . and display the results in fig. 6. 

Notice that, as expected, increasing decreases the error generation rate, while de- 
creasing Uk shortens the time to equilibration. Moreover, the scale of the error generation 
rate is correctly predicted by f lrlo]) . For instance, with u fc = 2.83^/, flrJo"]) sets an upper 
bound on the error generation near two percent per tp P , while the plot shows a maximum 
rate of approximately one percent per tp P . 



8. Discussion and Conclusions 

A hybrid algorithm for the accelerated simulation of Coulomb collisions has been 
presented. The algorithm is derived directly from the LFP equation without appealing 
to the ad hoc modeling used in other hybrid particle methods [ljj 24], permitting 
quantification of error sources and scalings. The accuracy and efficiency of the method 
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Figure 6: The systematic L dethermalization error as a function of time for three different 
initial conditions. 

is confirmed by the results of numerical simulations. Moreover, for this method the 
number of kinetic particles tends to zero as equilibrium is approached, thus recovering 
the efficiency of a fluid scheme for Maxwellian distributions, which is not a feature of the 
scheme in [6|. 

It is an unfortunate - but not unexpected - consequence of the approximations used 
to accelerate the algorithm that the hybrid scheme does not converge as At — > 0. It 
bears noting that one could make the hybrid scheme presented here convergent by tak- 
ing the relative entropy cutoff H c = 0(At k ) for any k > 2. However, such a scheme 
inevitably has the property that the speed-up factor S — > 1 as At — > 0, thus recovering 
the computational efficiency of a PIC scheme in the limit. The choice of k = 2 made 
in the preceding results is the unique choice that bounds the computational gain from 
below and the error from above. Moreover, we demonstrate in sections 6.3 and 7.2.4 the 
ability to predict the scale of the dominant error and find it to be small for the intended 
applications. 

For a given number of simulated particles, the entropy based scheme has been shown 
to be more accurate than the particle- and scattering angle-based schemes. However, the 
additional computational load incurred by the tracking of the passive scalars assigned to 
each particle is seen to effectively cancel this gain in the case of a spatially homogeneous 
relaxation process. However, in the application of the hybrid scheme to problems with 
other potentially destabilizing agents - e.g. spatial inhomogeneity or electromagnetic 
fields - we expect that the dethermalization error will become more significant, and it 
is reasonable to expect the entropy scheme's improved treatment of dethermalization to 
yield dividends. 

There are a number of directions in which the present work could be extended. The 
first is the treatment of spatially inhomogeneous problems, which is of obvious impor- 
tance for application to real world scenarios, and may reveal more benefits of the entropy 
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scheme, as just mentioned. A second is the extension to other collision operators, which 
only requires modification of the expressions for W and F. Thirdly, one might incor- 



porate unequal weighting of the simulated particles [18[, potentially including negative 
weights as in 17|. A fourth is to use the numerical evaluation of the expressions in ap- 
pendix C to evolve the relevant passive scalars, which may yield a more robust scheme. 
Other directions include the potential for adaptively choosing the number of moments 
used to approximate the particle distributions and the possibility of fusing multiple ki- 
netic particles when the relative entropy between them is small so as to further reduce 
computational cost, an approach similar to 15|, |24J. Each of these is a topic of current 
research of the authors. 
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Appendix A. Proof of Relative Entropy Theorem 

In 4.1, we state a theorem about decay of relative entropy: 

C(fJ m )\og(4-) dv<0. (A.l) 



Here, we provide a proof for the case of the Boltzmann collision operator, given by 

where da/dQ is the differential cross section of the inter-particle force mediating the 
collisions, the pre-collision velocities are v, v#, the post-collision velocities v', obey 
particle momentum and energy conservation, /' = f(v'), and similarly for the other 
evaluations of / and g. 

Because it adds little in the way of complication here, we treat the general case of 
species with distinct mass. This does require some additional notation: let the two 
species have distributions / and g, with Jm and gu being Maxwellian distributions for 
each respective species having identical mean velocity and temperature. 

We initially proceed as in the proof of the if -theorem. Let <p(v) be an arbitrary 
function of v. Writing out C(f,g), 

[ C{f,g M )ip{v)dv = [ B(0,\v-v*\)(fg' M *-fg M *)<p(-v)dndv*dv, (A.3) 

where B — |v — v*\da/dQ. Conservation implies 

m/v + m fl v* = m/v' + m g v*, (A. 4) 
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vrifV 2 + rrigvl = rrifV 12 + m g vf. (A. 5) 

We see that the operator in flA.3j) features one of the two symmetries exploited in proving 
the if-theorem: It is anti-symmetric with respect to exchange of the primed and un- 
primed velocities. Thus, we can write 

/ C(f,g M )<p(v) dv=- f B(9, |v - v*\)(f'g' M * ~ f9M*)(^P ~ V 3 ') dQdv*dv. (A.6) 

JR3 1 J 

Straightforward algebraic manipulation can be used to show that |v — v*| = |v' — v^|, even 
when m/ 7^ m g . A shorter argument can be made by appealing to the time symmetry of 
binary collisions. 

Now set <f = log(///jvf), so that 

<p-<ft = log(///') - \og(f M /f M ). (A.7) 

Now, note that a simple calculation and observation of (IA.4D and (|A.5[) implies that 

log(/ M ) + log(^) = \og{f' M ) + log(<4J. (A.8) 

This arises because each side is a linear combination of quantities that are invariant under 
collisions (total momentum, energy, and mass). As a result, we have log(/M / ' f' M ) = 
\og(g'M*l 9m*)- Plugging this into the formula for ip — <p', and then putting that in (IA.6I) . 
we have 

[ log (-£-) C(f, g M ) dv = l [ U'9'ku ~ f9M*)B log (^) dQdv*dv, (A.9) 

where we've now omitted the arguments on B. Now, since B is non-negative, it doesn't 
affect the sign of the integrand, and the rest is of the form {x — y)\og(y / x) , which is 
non-positive and zero just in case x = y. Thus, the entire integral is at most zero. 

It just remains to show when equality is achieved. As already noted, (x — y)\og(y/x) 
is zero just in case x = y, so the right side of flA.9[) vanishes if and only if fgM* = / '9m*- 
Taking the logarithm of both sides gives 

log(/) + log(^) = log(f) + log(<4J. (A.10) 

By (IA.8p . / = fu satisfies this equation. If we write h = f / fu-, (lA.lOj) reduces to h — hi . 
Fixing v, we can pick and v'^ such that v' has any value we like. Thus, the only 
solution to h — h! is h — c(x), independent of velocity. Since we assume (/) = (/a/), we 
must have c = 1. □ 



Appendix B. Results on Collisional Moments 



To compute the integrals in (|34p and (j3"5|) . we begin by adopting the notation of 14 
for the LFP collision operator. We write 

9 ffi„ r, % 



CM/,/m) = -— ^R-D._£j, (B.1) 
where 

H= 22» D^|^, (B . 2) 
m av m z avav 
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and 



7 



e 4 log A 

87T£§ 



(B-3) 



We've assumed all the particles under consideration have common mass m and charge 
e. All subsequent results can be straightforwardly generalized to the case where these 
quantities differ between species. 

Moreover, both here and in appendix C, we find it convenient to work in the rest 
frame of /m, so that um = and UjM = u,. 



Appendix B.l. Approximating Uj 
We define 

F jM 



vC F p(fj, f M )dv, 



(B.4) 



so that (d t Uj)M = FjM. Using ( IB.lj) . we may integrate by parts and use properties 
relating R and D (see 14[ for more detail) to find 



jM 



fjRdv. 



™ Jr3 

Then, using (IB. 2[) and (EJ), we can write an explicit expression for R. 

2jn M v 



R 



m v 



d erf(x) . 
x : erf (x) 

ax 



(B.5) 



(B.6) 



where um is the number density associated /m, and x = v/v t M- Under the assumption 
that fj is Maxwellian, the integral in (1B.5j) can now be evaluated numerically in general, 
and analytically in two important limits. 

Firstly, as mentioned earlier, near the beginning of a hybrid simulation, fj is well 
approximated by a 5-function at the particle's actual velocity, which now coincides with 
its mean velocity uj. Making this approximation, we easily find that 



jM 



m2v tM U jM 



u 



jM- 



d erf([7. 



jM, 



dx 



eii{U jM ) 



(B.7) 



where UjM = u jM/vtM- We'll refer to this expression for Fj M as F^ M , which is valid 
when v tj -C v tM (i.e. Tj < T M ). 

Secondly, at late times in the simulation, we expect that UjM *C v t j, so that we may 
approximate fj = / m (v; 1, Uj,Tj) by 



fj*f m (v;l,0,Tj) (l + ^p)- 



(B.8) 



Plugging this approximation into (IB.5P gives (see 1J] for more detail) 

1 



jM 



TjM 



- u jM, 



where 



TjM 



3v^Fm 2 ( 



Vtj + v. 



2 ^3/ 2 

tM 



r 



16 7n M 

This expression for Fjm will be called F™ M and is valid when Uju <^ v t j. 



(B.9) 



(B.10) 
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Appendix B.2. Approximating T 
We define 

W jM = \mf |v - u/'Crrifr /at) dv, (B.ll) 
so that (d t Tj)M = 2WjM/3. By expanding the squared term in W^m, we find 

W jM = \rn [ v 2 C FP (fj, f M ) dv - mUj ■ F jM . (B.12) 

Having already computed expressions for F, it suffices to compute what we'll call W, 
defined by 

W' jM = \m f v 2 C FP (fj, fu) dv. (B.13) 

Again, by plugging in the definition of Cpp into (IB.13P and integrating by parts - twice 
this time - then using the definitions of R and D, we find 



W' M = K 



^d erf(x) erf(x) 



dx x 



fjdv, (B.14) 



where 

(B.15) 

mv tM 

and x = v/vtM- Again, the integral in (IB. 14[) can be evaluated numerically for a general 
Maxwellian /, and analytically in two important limits. 

We again treat the case in which v t j v tM . As before, we approximate fj by 
5 3 (v — Uj). The integral for Wj M is now easily evaluated, and when combined with (1B.7j) 
and (1B.12|) gives 

WjM = n G -^l = 2 ^ erf f^) , (B.16) 

UjM mVtM Uj M 

In analogue to the previous subsection, this expression for WjM will be referred to as 
Wj M and is valid when v t j v t M- 

Finally, we consider the UjM *C vtj limit, just as we did when approximating FjM- As 
before, we approximate fj by a Taylor series expansion, this time keeping 

f jX f m{ y-l,0,T j )(l + 2 -^- U M (B.17) 

V U tj U tj J 

(we ignored the last term before because it gave no contribution to the previous integral). 
Using this expression for / in (IB.14j) and integrating gives 



W, 



jM 



TjM 



2 



mu jM 



(B.18) 



Again in analogue to the previous subsection, we will denote this expression for WjM by 
WjM and it is valid when Uj M v t j. 
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Appendix C. Numerical Integration of Collisional Moments 

Because of the length of the expressions in this appendix, we abbreviate our notation 
by dropping the subscripts j and M wherever possible, and again work in coordinates 
where Um vanishes. That is, we let fj — > f, v t j —> v t , ujm = % —> u, Fjm — > F, and 
W jM ->• W. 

In appendix B, we derived integral expressions for F and W, the collisional rates of 
change of the first and second moments of /, respectively: 



p = 1 



m Jus 



W 



fRdv 

d erf(x) erf(x) 



dx 



x 



fdv, 



(C.l) 
(C.2) 



where expressions for R and k are given in flB.6j) and flB.15j) . respectively, and x is as 
defined in appendix B. We then derived asymptotically valid analytic expressions for 
the case when / is Maxwellian with either v t <C v t M or u < w ( . In this appendix, we 
further demonstrate that the above expressions for F and W can expressed in terms 
of easily computable, one dimensional integrals depending on only two non-dimensional 
parameters and physical constants for a general Maxwellian f. 



Appendix C.l. Simplifying W 

Let us work in spherical coordinates with the z-axis aligned with u, the mean velocity 
of /. Then, a Maxwellian / may be written as 



/ 



exp 



u 



exp 



v 2 + u 2 



exp 



2uv cos tp 



(C.3) 



We observe that every term in (C.2) is spherically symmetric except for the right-most 
term in (10.3[) . The angular integration is thus easily performed, giving 



2tt pn 



exp 



2uv cos ip 



sin (pd(pd6 = ^ Vt sinh 



uv 



2uv 



(C.4) 



Next, we define T = v t Ai/v t , so that v/vt = Tx and u/v t = TU, with U as defined in 
appendix B. We can now rewrite 



w 



k r 



xG(x] u, r) 



deri(x) erf(x) 



dx 



x 



dx. 



where 



G{x\ U, T) = exp (-r 2 (x - U) 2 ) - exp (~r 2 (x + U) 2 ) 



(C.5) 



(C.6) 



The expression for W in (1C.5|) is now a one dimensional integral that is straightforward 
to evaluate numerically. 



Appendix C.2. Simplifying F 

We work in the same coordinates, and the expression (10.3|) for / still applies. However, 
the vector R contributes to the angular integrals. We write R = (v/v)R(x), where 



R(x) 



2777. 



M 



mV tM X 



d erf(x) . 

x ; erifxj 

dx 



(C.l) 
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The angular portion of the integral for F now reads 

r [ © exp ( 2 T sy ) sin * ^ d °- (c - 8) 

All the components of this vector vanish except for that aligned with u, because the 
integration against 9 yields zero in the other cases. The component along u may be 
evaluated, giving 

where 

^ = 2r 2 t/a;. (CIO) 

We then write F as 

F = 7^ (f) j^ xR{x)ew ( " r2 ( " 2 + ^ 2)) ( coshW ~ ^) ^ (an) 

Again, this integral is now easily evaluated numerically. 

The expressions (10.5j) and (IC.lip are used to generate the black curves in fig. 1 and 
may in theory be used to generate two dimensional (T, U) lookup tables which can then 
govern the evolution of the particle temperatures and mean velocities. 
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